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Abstract 

Parametric inference for spatial max-stable processes is difficult 
since the related likelihoods are unavailable. A composite likelihood 
approach based on the bivariate distribution of block maxima has been 
recently proposed in the literature. However modeling block maxima 
is a wasteful approach provided that other information is available. 
Moreover an approach based on block, typically annual, maxima is 
unable to take into account the fact that maxima occur or not simulta- 
neously. If time series of, say, daily data are available, then estimation 
procedures based on exceedances of a high threshold could mitigate 
such problems. In this paper we focus on two approaches for com- 
posing likelihoods based on pairs of exceedances. The first one comes 
from the tail approximation for bivariate distribution proposed by 



Ledford and Tawn (1996) when both pairs of observations exceed the 



fixed threshold. The second one uses the bivariate extension (Rootzen 



and Tajvidi , 2006 ) of the generalized Pareto distribution which allows 
to model exceedances when at least one of the components is over 
the threshold. The two approaches are compared through a simula- 
tion study according to different degrees of spatial dependency. Re- 
sults show that both the strength of the spatial dependencies and the 
threshold choice play a fundamental role in determining which is the 
best estimating procedure. 

AMS 2000 subject classification: primary 62G32, secondary 62M30. 
Keywords: Composite likelihood, Extremal dependence, Generalized Pareto 
distribution, Spatial statistics. 



1 



1 Introduction 



Extreme value theory for multivariate data is well established and there are a 
lot of results which can be used to characterize extreme values distributions 
(see, for instance, 



Resnick, 1987 Beirlant et al. 2004) 



From a statistical point of view (Coles, 2001) there are mainly two ap- 
proaches to fitting multivariate models: block maxima and threshold ex- 
ceedances. Block maxima is the most used approach: assuming a paramet- 
ric family of distributions, the model parameters are estimated on a block 
maxima sample. In practice, bivariate distributions of maxima have been fre- 
quently considered and there are few papers that deal with higher dimensions 



(Tawn, 1988, 1990 Coles and Tawn, 1991, 1994). 



Recently, spatial extreme problems have received an increasing attention 

fll996b; 



in the literature (see, among others, Coles (1993); Coles and Tawn 



Casson and Coles (1999 


); 


Buishand et al. 


( 


2008 


); 


Bel et al. 


( 


2008 


); 


Padoan 


et al. 


(2010 


) and for recent reviews 


Bacro and Gaetan 


( 


2012 


); 


Davison et al. 



(2012)). In particular max-stable processes (de Haan, 1984) take a prominent 
role in modeling the spatial dependence because they are a natural extension 
of the multivariate extreme value distributions for dealing with spatial data. 



However, in spite of the relative large number of instances (Smith, 1990 



Schlather, 2002 Kabluchko et al. 2009), parametric likelihood inference for 



such models is not easy because the related likelihood is unavailable or too 
cumbersome to compute. For this reason a methodology based on partial 



specification of the full likelihood, the composite likelihood (Lindsay, 1988 
Varin et al. 2011), has gained popularity. For spatial extreme data, the 



composite likelihood (Smith and Stephenson, 2009 Padoan et al. 2010) is 



built by taking product of pairwise distributions. Then parameter estimates 
are obtained by fitting the model to block bivariate maxima vectors. 

In this paper we are concerned with estimation procedures based on ex- 
ceedances of a large threshold. In doing this we expect that the efficiency 
of the estimates can be improved in all situations where the paucity of data 
prevents the traditional approach based on the block maxima. For instance, 
rainfall data are frequently collected every day but for few years. If annual 
maxima are considered, we have got estimates by using few repetitions. A 
larger data set can be obtained by dealing with daily exceedances of a large 
threshold, even if constraints of stationarity can reduce the temporal window 
to some months per years. 

Our aim is not suggesting new models for explaining the spatial behaviour 
of exceedances (see, for instance, Buishand et al. , 2008 Turkman et al. , 2010) 



but proposing and comparing two different approaches for composing likeli- 
hoods based on exceedances. The approaches rely on different specifications 
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of the bivariate distribution of the exceedances. 



The first one comes from Ledford and Tawn (1996). These authors have 
proposed a tail approximation for the bivariate distribution valid when both 
components exceed a large threshold. Since the approximation is valid for 
values simultaneously larger than the threshold, a censored version of this 
approximation is proposed for pairs of data having one or two components 
below the threshold. 

The second approach we have considered is based on a result for multi- 
variate exceedances due to Rootzen and Tajvidi (2006). This result allows 
us to model exceedances over a threshold using a multivariate extension of 
the generalized Pareto distribution when at least one of the components is 
over the threshold. The result is asymptotic and leads to approximating 
the conditional distribution given that at least one component is over the 
threshold. 

The paper is organized as follows. Spatial max-stable processes are pre- 
sented in Section [2] The two bivariate models for the exceedances are de- 
scribed in Section [3] and the estimation procedure we propose is detailed 
in section |4j Section [5] presents a simulation study where the two thresh- 
old exceedances models are compared, according different degrees of spatial 
extreme dependencies. Finally, some perspectives are discussed in Section |6l 



2 Spatial max-stable processes 



When we deal with extreme values, the traditional approach considers the 
maxima over a number of replications using the results from the extreme 
value theory. This theory embraces the univariate and the multivariate case 
(see, for instance, Beirlant et al. , 2004, for a recent account). In the case of 
spatial data maxima of observations over a fixed period are frequently col- 
lected on a finite subset of sites and the aim is to characterize the stochastic 
behavior on an uncountable set of unsampled sites. For this, multivariate 
models are unsuitable and their natural extensions are the max-stable pro- 
cesses Qde Haan[[l984| . 

Let {Z(s),s £ S} be a stochastic process on an index set S C M. d . We 
assume that Z k (-), 1 < k < n, are n independent copies of Z(-). The process 
Z(-) is a (spatial) max-stable process if the multivariate distribution of p > 1 



sites Z(si) 



Z(Sp), with Si G S and % 



,p, satisfies the so-called 



max-stability property, i.e. there exist sequences a n (s) > 0, b n (s), n > 1, 
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such that 
Pr 



Z k {sj) 



<z f l<j<p) =Pr(Z( Sj )<z j; l<j<p) 



max 

l<k<r 

(1) 

From the definition, any finite distribution of a max-stable process is an 
extreme value distribution. In the sequel we will denote the p dimensional 
finite distribution of a spatial max-stable process as G slt _ iS (zi, . . . ,z p ). 

If we suppose that we have n independent copies of a process Y(-) such 
that 

Y k (Sj) - n {Sj) 



lim Pr 

n— >oo 



max 

Kk<n 



a n (Sj) 



< Zj] l<j<P 



G 



81, 



\zu 



1 z p) 



(2) 



for suitable sequences a n (s) > 0, (3 n (s), we say that the distribution of 
Y(si), . . . , Y(sp), i.e. F SI s , is in the attraction domain of the extreme value 
distribution G S1) _ >S . This fact will be denoted by F S1) _ >S G T>(G Sl) ^ tSp ). 
There are several representations for max-stable processes, see for in- 



stance Smith (1990); Schlather (2002); Kabluchko et al. (2009); Lantuejoul 



et al. (2011). Just for illustration, in the sequel we will concentrate on the 



Gaussian extreme value process (Smith, 1990). 

Let II be a Poisson process of intensity x~ 2 dxds on (0, oo) x S. Then, 

Z(s) = max xf(s-u), (3) 

(ic,u)en 



defines a stationary max-stable process, with unit Frechet marginal distri- 
bution, i.e. Pr(Z(s) < z) = exp(— 1/z), z > 0, provided that /(•) is a 
non negative function such that J s f(s)ds = 1. Choosing / as the bivariate 
Gaussian density 



f( S ) = (27T)" 



-1/2 



exp 



1 



s'E" 



where E is a covariance matrix, leads to the Gaussian extreme value process. 



Smith (1990) gives an intuitive physical interpretation of the process <\3h in 



terms of rainfall-storms i.e. x and u are, respectively, the size and the center 
of a 'storm' that spreads out according the density /. 

For max-stable processes analytical expressions for multivariate distribu- 
tions are quite cumbersome and explicit formulas are known mainly for the 
bivariate case. Smith (1990) derived the bivariate distribution of ^ for two 
different sites and Sj, namely 



G SijSj (zi,Zj) = exp 



2 



log 



$ 



Z 3 



+ 



log 



Z 3 



(4) 
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where <£>(•) denotes the standard normal cumulative distribution function and 
an = y/(si — Sj)'T,~ 1 (si — Sj). In the case of a marginal Gumbel distribution, 



"13 
i.e. 



Pr(Z(Y 



J 3> 

<z) 



exp{— exp(— z)}, the bivariate distribution is given by 



G, 



exp 



-e 2! $ 



+ 



z 3 (|) 



+ 



More recently Genton et al. (2011) have derived multivariate expression 



for the Gaussian extreme value process but in practice only the three-variate 
case is simple to deal with. 

This fact has restricted the application of the parametric inference based 
on the full likelihood and motivated approaches based on the composite like- 



lihood (Lindsay, 1988). Roughly speaking (see section |4]for more details) the 
composite likelihood is an estimating function obtained combining marginal 
or conditional densities. 

However, these new advances still require a lot of temporal information, 
because we need to calculate the maxima over several period for providing 
enough sample data. This weakness can be overcome if we use all observa- 
tions that exceed some fixed (high) threshold. For this reason in the next 
section we will illustrate two different ways for modeling bivariate threshold 
exceedances and we will exploit these representations for building suitable 
composite likelihoods. 



3 Models for bivariate threshold exceedances 



Univariate extreme values results (Pickands, 1975 Davison and Smith, 1990) 
justify to model the conditional distribution of the excesses Y(s) — u through 
the Generalized Pareto (GP) distribution 

Pr (XM_Z^. < Z \Y( S ) > wj w 1 - (1 + £zr 1/f , z > 0, 

provided that F a (-) G T>{G). Here £ is a real shape parameter and a u is a 
positive scale parameter, which depends on the threshold u. 

The GP distribution provides an approximation of the tail of the distri- 
bution function in the sense that for a large enough value of u, we have 

Pr - U < z^j re 1 - C {1 + izY 1 ^ , z>0, (5) 

where ( = 1 — F s (u). 
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Now we consider the bivariate random variable Y(si),Y(sj) with the same 
marginal distributions and F Si>s . 6 T>(G s . iSj ). According (J5|, the integral 
transformation 



y k = - io g a-c 



y(sjfc) - u 



(6) 



has unit Frechet distribution. 

The bivariate distribution F s . s .{yi,yj) can be approximated on (u, oo) 2 
(Ledford and Tawn, 1996) as 



F St , Sj (yi,yj) w G s . iS . 



(7) 



Formula ([7]) gives us an approximation of the bivariate distribution in the 
tail when both components are greater than threshold u. 

Another approach (Beirlant et al. , 2004 Rootzen and Tajvidi 2006) con- 
siders a bivariate distribution of large values, given that at least one of the 
components is large. This distribution is a bivariate analogue for the GP 
distribution. 



Provided that F. 



■) G V(G Si>Sj (; •)) and < G< (0,0) < 1, Rootzen 



and Tajvidi (2006) show that for large u 

Y( S< 



Pr 



Y s, 



u 



< Zi 



u 



< Zj I Y(si) > u or Y(sj) >«) ~ H s ^ Sj (z i} Zj) 



where 



/ / (.Tj . Xj 



log 



log(G Si)S .(0, 0)) G SuS . (min(a; i , 0), mm.(xj, 0)) ' 



Actually the Rootzen-Tajvidi's result applies to a dimension greater than 
two, and in the sequel we refer to the distribution (J8|) as the bivariate GP 
distribution. 



4 Inference procedures for spatial high thresh- 
old excesses 

As mentioned in section [2j likelihood inference for spatial max-stable pro- 
cesses is difficult because we don't known an analytical expression for the 
multivariate distribution of the majority of the spatial max-stable processes. 



The composite likelihood (CL) (Lindsay, 1988) is an inference function 
built by multiplying the likelihood of marginal or conditional events. More 
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precisely, let {/(■, 9), 9 G B} be a parametric family of joint densities for the 
observations yi,...,y n G D C W 1 and consider a set of events {Ai : C 
e J C N}, where 5 is a a-algebra on D. The logarithm of the CL is 
defined as ( Lindsay , |1988[ ) : 



d (°) = E Wi lo S/((2/i» • • • ) 2/n)' e ^ 



where Wi are non negative weights to be fixed. 



For max-stable processes, Padoan et al. (2010) have proposed a CL ap 



proach based on the pairwise distributions for estimating the parameters of 
the distribution of random vectors of block maxima. Their estimation pro- 



cedure has been implemented in the R package Spat ialExtr ernes (Ribatet 



2011) available on the CRAN repositories. Here we borrow the same idea 



but we use exceedances instead of block maxima. 

We assume that data {y t \, . . . , yt P }, t = 1, . . . , T, are T independent real- 
izations of the vector {Y(si), . . . , Y(s p )}, where si, . . . , s p are p sites belong- 
ing to S. 

Let f Si , Sj (-, - ]9) be the bivariate density associated to the distribution ff 
or <\8n. The CL estimate 9t is obtained maximizing 



T p p T 

mo) = E E E w » ] °g f«* (w*' v& e ) = E u *w 
t=i t=i j>i 



(9) 



where Wij > designate user-specified weights (Padoan et al. 2010) and 
Ut{9) = y>; ,yj;. l a;,U^.U^(Y l ,.) u :9). 

Note that there is a fundamental difference between the estimating func- 
tions, when we choose (|7j) or ([8]). In using ^ we exploit a tail approximation, 
moreover the approximation is sound on the region (u, oo) 2 and does not ap- 
ply to pairs of observations outside that region. For pairs such that at least 
one component does not exceed u, we use the censored approach as proposed 
in 



Ledford and Tawn (1996). The likelihood contribution of each pair is 

V o6 G s . iSj (y ti , y tj ; 9) if y ti > u, y tj > u 
(y u ,u;9) i£y ti >u,y tj <u 



fsi,sj {yti-, ytj', 9) 



V a G s 



G SuS .(u,u\ i 



if y ti < u, y t j > u 
if y ti < u, y tj < u 



where 9 = (a u ,£,fl) and is the (possibly) vector of spatial dependence 
parameters. Instead of maximizing cl{9) with respect to all parameters, we 
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can estimate 9 in two steps. First of all we estimate the marginal parameters 
£ et a u , then their estimates, | and d u , are plugged in ^ for obtaining an 
estimate for /3 maximizing cZ(£, a u , 0). 

When we use ([8]), differently we keep the pairs such that there is at least 
one exceedance and the number of pairs in ^ is a random variable. Moreover 
we refer to conditional events. The associated density is given by 



fsuSjiVtij Vtj'-i 9) 



1 



log G SijSj (u,u; 9) 
V a bG Si , s .(y ti ,y tj ;9) 



V a G s<jS Ayu, Vtj] 9)V b G Si , s (y ti , y tj ; 6)' 



G s . tS .(y t i,y tj ;9y 



for (y u , ytj) 4- (—oo,u] 2 . In this case the components of the parameter vector 
9 are the marginal parameter a u and the dependence parameters /3. 

Exceedances above u will appear in spatial clusters and in this respect 
the weights Wij in ^ should handle pairwise dependencies inside and outside 
of the clusters. 

Composite log-likelihood, as linear combination of log-likelihoods, allows 
to obtain unbiased estimating equations under classical regularity conditions 



(Heyde, 1997). So, in principle, it would be possible to come up with an 



optimal way of combining the individual scores V^log/ Si)S .(y^, ytf, 9) and 
determining optimal weights (Kuk, 2007). In the next section , more sim- 



ply, we look at weights for reducing the computational burden as in |Padoan 
et al. (2010), like as the cut-off weights, namely u>„- = 1 if 



\Si S n 



<8, 

and otherwise, for a fixed 5 > 0. There is some evidence for Gaussian 



random processes (Davis and Yau, 2011 Bevilacqua et al. 2012) that this 



choice improves not only the computational efficiency but also the statistical 
efficiency 



Under suitable conditions (see the Appendix in Padoan et al. , 2010) the 
maximum composite likelihood estimator for 9 can be proved consistent and 
asymptotically Gaussian, for large T. The asymptotic variance is given by 
the inverse of the Godambe information matrix 



g T (9) = n T (9)[M9)]- 1 n T (9) 



(10) 



where H T {9) = E(-V 2 cZ T (fl)) and J T {9) = Var(Vcl T {9)). 

Standard error evaluation requires consistent estimation of the matrices 
Ht and Jt- Assuming strong-stationarity in time, we get 



U T {9) = -TE{V 2 f/i(#)} 
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that can be estimated by Ht = — V 2 c/t(^t)- Estimation of the matrix 

T T 



t=l k>t 



requires some care. When we deal with real data, for instance environmental 
data, exceedances are seldom independent in time. In such case declustering 
techniques (Nadarajah, 2001) could be used to overcome the dependence in 



time. An alternative way is employing a subsampling technique ([Carlstein 

ti 
1 



1986]) to estimate Jt- The subsampling method consists of estimating Jt 

,M, of 



over M overlapping temporal windows Dj C {1, . . . , T}, j 
size dj by using the expression 



Ji 



T A Vcl Dj 



)VcZ 



D, 



7T 



Al 



3=1 



d J 



where Vein is the composite likelihood score evaluated over the window Dj. 
An estimate of the asymptotic variance is then given by 



Hrp ^ Jj'Hrp ^ . 



5 Numerical results 



In this section we describe a simulation study with the aim of examining and 
comparing the performances of CL estimates for the parameters of spatial 
max-stable processes. Three versions of the CL are contrasted: two based on 
threshold exceedances that use the bivariate distributions ([TJ) and (|8~]), noted 



LT and RT, respectively, and one based on block maxima data (Padoan et al. 



2010), PRS hereafter. The advantage of threshold methods over the block 



maxima one is well known even if this advantage deserves to be evaluated 



(Madsen et al. 1997). 



set S 



We have considered the G aussian extrem e value process (|3|) and we have 
200 



150 



150 
300 



as m 



Padoan et al. 



(2010). The spatial process is 



observed on n 2 points located over a n x n regular grid {k, . . . , n* k} 2 , where 
n = 5,7. Here we present only the results for the 49 sites, because evidences 
are virtually the same with the smallest grid (25 sites). Three lags k for the 
spatial sites, namely k = 1,5, 10, are considered in order to take into account 
different levels of the spatial extremal dependence. 
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The spatial extremal dependence of the maxima can be characterized 



using the extremal coefficient function (Smith, 1990; Schlather and Tawn 



2003|), v(-), given by 

Pr (max(Z(s), Z(s + h)) < z) = exp{-v{h) / z). 

For any vector h, we have 1 < v(h) < 2 and v (h) = 1 corresponds to 
the perfect dependence, instead for independent maxima we have v(h) = 2. 
In the case of the Gaussian extreme value process, the extremal coefficient 
function is 

v(h) = 2$ I — 

Figure [T] shows the values of this function for the grids over 49 spatial sites and 
for the different spatial lags. As expected, the spatial extremal dependence 
decreases as the spatial lag increases. For k = 1, the extremal dependence 
keeps strong on the whole grid. For k = 10, the dependence appears strong 
only for neighbouring sites. 

In our experiments, for each simulation, we have considered 40 years and 
for each year we have a sample of 91 independent daily observations. Our 
aim is to mimic a framework where maxima are taken over a fixed season 
of the year. Then we have considered 500 simulations for evaluating mean, 
standard deviation, and root mean square error of the estimates. 

Simulations have been carried out using the rmaxstab function of the R 



package SpatialExtremes (Ribatet, 2011). This function allows to simulate 



a finite realization (Z(s\) } . . . , Z(s p ))' of a Gaussian extreme value process 
with unit Frechet marginal distributions. In order to fulfill the condition 



< G SijSj . (0, 0) < 1 of Theorem 2.1 in Rootzen and Tajvidi (2006), we 
have transformed the data to the Gumbel scale, namely Z*(st) = \og(Z(si)), 

1 = 1, . . . , n. The threshold u for all sites has been set equal to the empirical 
0.98-quantile, calculated over all data in each Monte Carlo simulation. 

As mentioned at the end of Section |4j a careful choice of the weig hts in Q 
could improve the statistical efficiency in addition to the computational one. 
In fitting the models, we have considered three sets of weights, according to 
different values of S, where 5 has been set equal to the a-quantile, q a , of the 
distances among all pairs of sites and a = 0.25, 0.50, 1.00. 

Figures ([2]), ^ and ^ contrast the variability of the estimates for grids 
with different spatial lags and different set of weights derived from 5 = go.25 
and 5 = gi.oo- The RT estimates always display the smallest variability but 
raise some concerns according different degrees of spatial dependence. When 
we use all the pairs (5 = gi.oo) m forming the CL function, the bias of the 
an and cr 2 2 estimates seriously increases in case of weak spatial dependency 
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(spatial lag k = 10). On the other hand LT results point out an efficiency 
similar to the PRS method and the bias of the estimates keep within the 
reasonable bounds, regardless of the spatial dependence and the number of 
the pairs. 

In Tables [T], [2] and [3] we measure the efficiency of the estimates in terms of 
the root mean square errors (RMSEs) under different extremal dependence 
and different amount of pairs. Bold figures highlight the minimum value of 
RMSEs. When the extremal dependence is strong (grids with lag k = 1), 
the RT method clearly outperforms the LT one for whichever value of S. For 
moderate spatial dependencies (grids with lags k — 5), there is again an 
advantage in prefering the RT method with respect to the LT. Results for 
weak spatial dependencies (grid with lag k = 10) are quite differents and 
put forwards the role of the number of pairs: for 5 = qo.25, RT approach 
gives smallest RMSE than LT and for 5 = gi.ocb the converse happens. For 
8 = <?o.50, there is not clear indication of which method behaves better. 

However our results, consistent to the literature on weighted composite 



likelihood (Joe and Lee, 2009 Davis and Yau 2011; Bevilacqua et al. 2012), 



display that including many pairs in the pairwise likelihood can harm the 
CL estimator, suggesting that we should retain as few pairs as possible. 

We have already remarked that the RT estimates always had smaller 
standard errors than the LT ones. Nevertheless, they also lead to a larger 
RMSEs when the spatial dependence is weak and we use all the pairs because 
a large bias is present. This observed bias in the RT estimates seems also 
depend on the choice of the threshold. As matter of proof, we did the same 
Monte Carlo experiment setting the threshold u equal to the empirical 0.90, 
0.95 and 0.98 quantiles of the data. These values correspond to the usual 
choice in the real data analysis. Table [4] reports the bias of the parameter 
estimates according the three values of the threshold. We can see that for 
lower thresholds, the absolute value of the bias increases for the an and 022 
estimates following the RT approach. On the other hand the bias of the LT 
estimates is always moderate. 



6 Discussion 



In this paper, we have compared two threshold exceedances models to esti- 
mate the parameters of a particular instance of spatial max-stable processes, 



namely the Gaussian extreme value process (Smith, 1990 Padoan et al. 



2010) of spatial max-stable processes. Each model is based on a particular 



bivariate extreme value property (Rootzen and Tajvidi, 2006 Ledford and 



Tawn, 1996). However our estimation strategy can be easily extended to 
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other spatial max-stable processes (Schlather, 2002 Kabluchko et al. 2009), 



provided that the expression of the bivariate distribution is easy to evaluate. 

Our simulation experiments point out that both methods are valuable and 
choosing between the RT and LT approach relies on the degree of the spatial 
dependence. When the spatial dependence between high values is strong, the 
RT approach seems preferable. The RT approach for weak spatial dependent 
data displays a significative bias in the estimation. On the other hand LT 
approach does not suffer from this problem. 

A possible explanation of this evidence is that the bivariate distribution 
8b deals with pairs of observations with at least one component exceeding the 



threshold u. As Rootzen and Tajvidi (2006) have underlined in their paper, 
the multivariate generalized Pareto is degenerate in the case of independence, 
so it is not surprising that the weak extremal dependence can spoil the CL 
estimates. In this regard, our findings are consistent with the results obtained 
by Rakonczai and Tajvidi (2010) for a bivariate logistic extreme value with 

is not the true conditional 



weak dependence. Moreover, the distribution 
bivariate distribution but only its asymptotic approximation. Again in the 
case of weak extremal dependence, the convergence may be slow, leading to 
worry regarding to the usual behaviour of the CL estimates. 

Nevertheless our simulations highlight that a clever choice for the weights 
of the CL could significantly improve the efficiency of the estimates. This is 
especially true under weak spatial dependency, where the RT approach gives 
the best RMSE results. Indeed, when we consider all the pairs, pairs of sites 
which are not spatially dependent could be linked because of an exceedance 
of one component. In the bivariate CL framework, this should lead to a 
misleading link in the extremal association and, as a consequence, contribute 
to an incorrect estimate. 

This fact entails that there is a effective need of simple practical rules 
for fixing the weights. However rules based only on the extremal coefficient 
function appear as not really adapted and more theoretical work seems nec- 
essary. A possible solution that we will explore in the future is to choose the 
weights by minimizing a certain norm of the Godambe information matrix 
(10) as in Bevilacqua et al. (2012). 
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k = 5 k = 10 




Figure 1: On the left top corner: an example of a simulation from a Gaussian 
extreme value process, with a n = 200, a 2 2 = 300, a 12 = 300. Three different 
regular grids, with different spatial lags {k = 1,5,10), are superimposed. 
The other figures are the plots of the extremal coefficient function v(h). 
This function is evaluated at the distances among 49 sites on regular grids 
with different spatial lags. Solid lines (dashed lines) correspond to the value 
8 = %.25 (8 = 9o.50, respectively). 
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Parameter 


o-ii 


022 


012 




True value 


200.00 


300.00 


150.00 


5 = go.25 












RT 


209.13 (29.62) 
30.99 


308.71 (44.14) 
44.99 


156.49 (30.06) 
30.75 




LT 


204.24 (54.85) 
55.01 


301.80 (80.63) 
80.65 


152.71 (46.81) 
46.89 




PRS 


216.94 (47.77) 
50.69 


320.01 (67.19) 
70.10 


162.03 (46.94) 
48.46 


5 = 90.50 












RT 


209.00 (29.18) 
30.54 


308.46 (43.53) 
44.34 


156.31 (29.57) 
30.24 




LT 


204.75 (54.35) 
54.56 


302.42 (79.73) 
79.76 


153.08 (46.53) 
46.64 




PRS 


217.53 (49.03) 
52.07 


320.57 (68.79) 
71.80 


162.36 (48.16) 
49.73 


5 = gi.oo 












RT 


208.79 (28.09) 
29.44 


308.29 (42.12) 
42.92 


156.02 (28.45) 
29.07 




LT 


204.88 (54.63) 
54.85 


302.43 (79.90) 
79.93 


153.13 (46.79) 
46.89 




PRS 


218.66 ( 51.33) 
54.61 


321.70 (71.49) 
74.71 


163.00 (50.13) 
51.79 



Table 1: Simulation results for a grid of 49 sites with spatial lag k — 1. For 
each estimation approach (RT, LT, PRS), the mean and the standard errors 
(in brackets) of the estimates are reported on the first line. RMSE values are 
given on the second line. 
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Parameter 


on 


022 


012 




True value 


200.00 


300.00 


150.00 


5 = go.25 












RT 


204.41 (16.33) 
16.91 


305.37 (25.63) 
26.19 


152.06 (16.92) 
17.04 




LT 


201.06 (35.39) 
35.41 


302.03 (56.71) 
56.74 


151.55 (32.59) 
32.63 




PRS 


206.78 (29.79) 
30.55 


312.02 (49.35) 
50.79 


155.90 (31.49) 
32.04 


<5 = <?0.50 












RT 


205.58 (15.54) 
16.51 


304.47 (24.01) 
24.42 


148.33 (15.89) 
15.97 




LT 


201.53 (35.92) 
35.95 


302.52 (57.17) 
57.23 


151.87 (33.07) 
33.13 




PRS 


208.22 (33.67) 
34.66 


313.99 (54.34) 
56.11 


156.79 (35.34) 
35.99 


5 = 9i.oo 












RT 


216.32 (13.80) 
21.37 


306.31 (19.81) 
20.79 


133.11 (13.26) 
21.47 




LT 


202.20 (37.21) 
37.27 


303.47 (58.34) 
58.44 


152.31 (33.65) 
33.73 




PRS 


211.34 (41.89) 
43.40 


318.21 (64.93) 
67.43 


158.78 (42.98) 
43.87 



Table 2: Simulation results for a grid of 49 sites with spatial lag k = 5. For 
each estimation approach (RT, LT, PRS), the mean and the standard errors 
(in brackets) of the estimates are reported on the first line. RMSE values are 
given on the second line. 
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Parameter 


on 


022 


012 




True value 


200.00 


300.00 


150.00 


5 = go.25 












RT 


207.38 (7.24) 
10.33 


299.57 (11.98) 
11.98 


138.27 (7.98) 
14.18 




LT 


199.91 (18.79) 
18.77 


300.53 (31.43) 
31.41 


150.21 (18.07) 
18.07 




PRS 


201.55 (19.56) 
19.60 


305.64 (32.64) 
33.09 


152.30 (20.97) 
21.08 


<5 = <?0.50 












RT 


225.64 (6.59) 
26.47 


309.70 (10.52) 
14.30 


122.79 (7.06) 
28.11 




LT 


200.16 (19.56) 
19.56 


300.98 (32.61) 
32.61 


150.39 (18.65) 
18.65 




PRS 


202.43 (23.68) 
23.78 


307.70 (39.42) 
40.12 


153.23 (25.76) 
25.93 


5 = 9i.oo 












RT 


315.27 (7.45) 
115.51 


397.33 (10.25) 
97.87 


120.22 (7.19) 
30.63 




LT 


200.56 (22.75) 
22.75 


302.21 (37.72) 
37.75 


150.96 (21.70) 
21.70 




PRS 


204.09 (31.00) 
31.24 


311.48 (52.11) 
53.31 


155.00 (33.83) 
34.17 



Table 3: Simulation results for a grid of 49 sites with spatial lag k = 10. For 
each estimation approach (RT, LT, PRS), the mean and the standard errors 
(in brackets) of the estimates are reported on the first line. RMSE values are 
given on the second line. 
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Table 4: Bias values for the RT and LT estimates using all pairs and three 
different threshold values, namely set to the 0.9, 0.95 and 0.98 empirical 
quantile. The spatial lags k are taken equal to k — 5 and k = 10. 



k = 5 k = 10 





0.9 0.95 0.98 


0.9 0.95 0.98 


a u RT 


50.48 30.65 16.32 


226.47 167.81 115.27 


LT 


0.80 1.55 2.20 


0.06 -0.04 0.56 


(T22 RT 


40.03 19.27 6.31 


222.36 155.91 97.33 


LT 


0.99 1.89 3.47 


0.18 0.26 2.21 


(7i2 RT 


-16.48 -18.61 -16.9 


-8.3 -20.17 -29.78 


LT 


0.65 1.33 2.31 


0.10 0.07 0.96 
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5 = <?0.25 



5 = 91.00 



S - 




RT LT PRS RT LT PRS 



Figure 2: Box-plots of the estimates of an- The true value for an is 
indicated by a dotted line. Each row in the panel indicates the results over 
a grid with lag k — 1,5, 10, respectively. In the first (resp. second) column 
of the panel the composite likelihood weights correspond to 5 = g .25 (resp. 
$ = Qi.oo)- 
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5 = <?0.25 8 = 91.00 




Sj - 



Figure 3: Box-plots of the estimates of 022 ■ The true value for 022 is indicated 
by a dotted line. Each row in the panel indicates the results over a grid with 
lag k — 1,5, 10, respectively. In the first (resp. second) column of the panel 
the composite likelihood weights correspond to 5 — g .25 (resp. 5 = 91.00) • 
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Figure 4: Box-plots of the estimates of au- The true value for au is indicated 
by a dotted line. Each row in the panel indicates the results over a grid with 
lag k — 1,5, 10, respectively. In the first (resp. second) column of the panel 
the composite likelihood weights correspond to 5 = g .25 (resp. 5 = gi.oo)- 
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